In this document, I want to provide a write-up of how this model is similar to and different from the versions previously implemented. This model is meant to provide assembly of multiple sites at the same time, that may, or may not, be connected in some fashion. Along the way, I will mention some of the input and output format that I am expect as documentation. At the end, I will be presenting some preliminary results. I especially want to use this as a vehicle to think about how to analyse those results.
library(dplyr) # Data Manipulation
library(ggplot2) # 2-D Plot
library(plotly) # 3-D Plot
library(ggfortify) # used for biplots of PCAs
library(vegan) # Ecological analysis mega-package
# library(RMTRCode2)
# https://stackoverflow.com/a/7172832
ifrm <- function(obj, env = globalenv()) {
obj <- deparse(substitute(obj))
if(exists(obj, envir = env)) {
rm(list = obj, envir = env)
}
}
First, we load up the preliminary results. In this case, we are loading a system in which 34 basal species and 66 consumer species form the pool for 10 unconnected environments. The pool and interaction matrix were assembled with the default parameters from Law and Morton’s 1996 work.
load(file.path(
"..", "experiments", "MNA-FirstAttempt-Result-Env10-None.RData")
)
In total, 9320 events were used in these environments, with the species and environment invasion both randomly assigned. The number of arrival and extinction events were controlled to both be half of this number. We chose this number due to the coupon collecting problem. In particular, we use the result that the probability of encountering each species is bounded: \[\text{Pr}(\text{Draws} < n \log_{e} n + c n) \rightarrow \exp(-\exp(-c)) \text{ as } n \rightarrow \infty\] where \(n\) is the number of species and \(c\) is a constant. For our purposes, we choose \(c = 5\) so that we have a probability of about \(99.3\%\) of seeing each species in each environment. In practice, we failed to observe 2 species-environment combinations. Notably, nearly every species had at least one successful invasion; 72 did not.
The initial abundance was set to be 4000 times the elimination threshold, in line with work on minimum viable populations (Traill et al. 2007). The elimination threshold is admittedly more arbitrary, since it sets an effective individual-area relationship. For this calculation, I set it to \(10^{-4}\), in line with our previous calculations. This is large enough to avoid numerical difficulties from precision, while being low enough to represent a decent sized region.
Since each population is assembling simultaneously, I chose to use exponential waiting times for the inter-species arrival and extinction times. Note that these rates are shared between species and environments, but arrival and extinction are fully independent of each other. Species and environment affected in each event was chosen uniformly at random. The question then is how to set the rate.
To set the rate in this case, I chose to set it to the largest eigenvalue magnitude of the per-environment interaction matrices. This magnitude corresponds to the strongest response we can see from the interaction matrix and, if the interaction matrix is a good approximation for the Jacobian around a stable fixed point (which is not guaranteed), indicates the characteristic time scale of the decay to equilibrium. Hence, (overall) arrival rates and (overall) extinction rates should happen on the same timescale as the (largest) dynamics in the system. Since there are 10 environments, we should then expect that 10 characteristic time scales, on average, should occur in between arrival events in the same environment.
With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.
LawMorton1996_PlotAbundance(result$Abundance[, c(1, 2:101)]) -> obj;
obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)
Every vertical line is a species introduction or extinction by neutral dynamics.
Perhaps more intriguing might be some sense of the biodiversity that we have in each system. We break the abundance results up by environment, then calculate the number of non-zero abundance curves at each time point. We also calculate the Shannon entropy (reminder: higher entropy means more uncertainty which means a flatter distribution).
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(ifelse(x != 0, x * log(x), 0))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Species = species,
Environment = i)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
Evenness = Entropy / log(Richness)
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Richness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Richness = mean(Richness)
),
mapping = ggplot2::aes(
x = Time,
y = Richness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
So richness hovers around a similar regime throughout the majority of the simulation. Note in this plot that we have emphasised one environmental curve and superimposed the mean in black. We do manage to reach heights of 11 species in one environment, but these heights are shortlived. Instead, we seem to observe a (time and environment averaged) value of 8.3677307. (If we consider the first 10,000 time units as burn-in, we instead see a value of 8.2536852.)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Entropy,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Entropy = mean(Entropy)
),
mapping = ggplot2::aes(
x = Time,
y = Entropy
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Entropy has a similar, but highly erratic, behaviour. If one tries to follow one of the entropy curves, then one sees that they have fairly substantial periods of almost smooth behaviour followed by suddenly very noisy behaviour, and noise seems to be the dominant mode if one tries to examine the low alpha environment in the background. There are some easy to make predictions. Extinctions reduce the entropy in the system, as you become more certain about what remains. Analogously, arrivals increase the entropy. We can probably better see the relationship beyond these principles by plotting entropy against richness and connecting observations by environment and time.
# ggplot2::ggplot(
# Diversity %>% dplyr::filter(Environment < 4),
# ggplot2::aes(
# x = Richness,
# y = Entropy,
# group = Environment,
# color = Time
# )
# ) + ggplot2::geom_path(
# ) + ggplot2::guides(
# alpha = "none"
# )
plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
mode = "lines", opacity = 1, line = list(color = ~Time))
It seems quite hard to tell, but there does not appear to be any particular orientation (clockwise, counter-clockwise) or similar pattern here.
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Evenness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Evenness = mean(Evenness)
),
mapping = ggplot2::aes(
x = Time,
y = Evenness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Evenness helps highlight that this is a high variance process but with a relatively constrained mean.
We can, of course, flip the idea on its head. Instead of examining the diversity of species within environments, we can look at the diversity of environments occupied by species. Since very few species end up occupying environments, we just look at richness. Unfortunately, this is quite a memory exhaustive task.
EnvDiversity <- lapply(
1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
environments <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Abundance = abundSum,
Species = i,
Environments = environments)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
EnvDiversity <- dplyr::bind_rows(EnvDiversity)
ggplot2::ggplot(
EnvDiversity %>% dplyr::filter(Richness > 1),
aes(x = Time, y = Richness, color = Species)
) + geom_point(
alpha = 0.01, size = 3
) + guides(
color = "none"
)
One immediately interesting trend here is that very few species are present across more than 5 environments at a given time. Indeed, only 2, 4, 5, 12, 13, 14, 22, 24, 29, 30, 32, 33, 35, 37, 43, 44, 47, 55, 57, 61, 68, 69, 71, 72, 76, 77, 78, 86 are ever present in more than 5 environments at once. We can also examine how long these periods occur for by species by tabulation. We block times so that entries are all of the same unit length, and round the average richness during the given time unit.
with(EnvDiversity %>% mutate(
Time = floor(Time)
) %>% group_by(
Time, Species
) %>% summarise(
Richness = round(mean(Richness)),
.groups = "drop"
),
table(Species, Richness))
(For the desktop, the amount of data we generated is too much, so we need to reduce the amount of data we use. To do so we average over time blocks, here of length 100.)
AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
We can perform a PCA and see if are data can be summarised by a small number of dimensions. As we shall see, constraining to the first 25 principal components does not harm the system much.
PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
There is not actually a lot of dependence within the system it appears. Consider the amount of variance explained by the first six principal components (ordered, as is tradition, by amount of variation explained).
head(summary(PCA)$importance[2, ])
The amount of explained variation is as follows.
sum(summary(PCA)$importance[2, ])
The traditional biplot follows, but despite the seeming presence of patterns, so little of the variance is explained that is probably not worth further examination.
ggplot2::autoplot(PCA, loadings = TRUE,
data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
We can try again with presence-absence data instead.
AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::mutate(
dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCAPA)$importance[2, ])
So it is a bit better, but not by much.
ggplot2::autoplot(PCAPA, loadings = TRUE,
data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
I should note that this is not a failure of the method persay. What this says to me is that dimension reduction of the system cannot reduce the system to a human-readable set of descriptors. The system is still being reduced (from 100 Species * 10 Environments to 25 Components to cover
sum(summary(PCA)$importance[2, ])
of the variation). My initial hypothesis for when the systems are coupled is that, as coupling increases, the number of principal components should decrease, since one system’s changes will be better predictors of the next system’s changes. (An interesting related question; do the number of principal components correlate with the amount of biodiversity in the system?
ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)
Since trying to reduce the system dimensionally does not work (well enough) a next attempt might be to consider how the pair-wise beta diversity changes over the course of the simulation. Initial problems include that there are quite a few measures of beta diversity as well as that the (square of the) number of environments will determine how many entries we need to consider.
To try this, we are going to reorganise our data into data frame with attributes of environment, time, and traits/species.
Environments <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
retval <- data.frame(Environment = toString(i),
Time = time
)
retval <- cbind(retval, env)
colnames(retval) <- c("Environment", "Time",
paste0("Basal", 1:34),
paste0("Consumer", 35:100)
)
return(retval)
},
abund = result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
),
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)
One perhaps interesting idea is to try to group the environments by cluster by considering the abundances (or presence-absence) of the species present as traits. The question then is whether the environments have a tendency to attract to specific points or if they wander around each other without relation. (The latter is more neutral.) If they appear to be convergent (which so far would seem to disagree mostly with the alpha diversity analyses), then that would imply that dynamics determine the majority of the system’s fate, while if they instead wander more randomly, then they would appear to be dominated by the neutral mechanisms. (Of course, this is affected by the rate of the neutral mechanisms, so it exists on a definite continuum.) (Note, of course, that cluster analysis usually includes something like PCA, so we would not necessarily expect a different result.)
My working hypothesis for how to use the distance measures is that we need a metric measure (Kindt and Coe, 2005) with any identification scrubbed off (site, time collected). Kindt and Coe suggest that taking the square root of Bray-Curtis makes it metric and thus usable for mathematical distance analyses (such as hierarchical clustering). The vegan package recomends using Jaccard instead of Bray-Curtis, for the same reason.
EnvironmentDistance <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard")
EnvDistClust <- hclust(EnvironmentDistance)
plot(EnvDistClust, labels = FALSE)
Note that are major breakdowns extremely early yielding a large number of clusters but there are expected to be 10 environments (if historical contingency mattered most) or a small number of clusters (if they are converging). The large number of clusters says to me that the neutral dynamics are jumping rapidly from cluster to cluster.
We can try again with presence absence instead.
EnvironmentDistancePA <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)
EnvDistClustPA <- hclust(EnvironmentDistancePA)
plot(EnvDistClustPA, labels = FALSE)
We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.
load(file.path(
"..", "experiments", "MNA-FirstAttempt-Result-Env10-Line.RData")
)
With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.
LawMorton1996_PlotAbundance(result$Abundance[, c(1, 2:101)]) -> obj;
obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)
Every vertical line is a species introduction or extinction by neutral dynamics.
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(ifelse(x != 0, x * log(x), 0))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Species = species,
Environment = i)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
Evenness = Entropy / log(Richness)
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Richness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Richness = mean(Richness)
),
mapping = ggplot2::aes(
x = Time,
y = Richness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Entropy,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Entropy = mean(Entropy)
),
mapping = ggplot2::aes(
x = Time,
y = Entropy
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 1109 row(s) containing missing values (geom_path).
Warning: Removed 104 row(s) containing missing values (geom_path).
plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
mode = "lines", opacity = 1, line = list(color = ~Time))
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Evenness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Evenness = mean(Evenness)
),
mapping = ggplot2::aes(
x = Time,
y = Evenness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 1440 row(s) containing missing values (geom_path).
Warning: Removed 136 row(s) containing missing values (geom_path).
EnvDiversity <- lapply(
1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
environments <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Abundance = abundSum,
Species = i,
Environments = environments)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
EnvDiversity <- dplyr::bind_rows(EnvDiversity)
ggplot2::ggplot(
EnvDiversity %>% dplyr::filter(Richness > 1),
aes(x = Time, y = Richness, color = Species)
) + geom_point(
alpha = 0.01, size = 3
) + guides(
color = "none"
)
with(EnvDiversity %>% mutate(
Time = floor(Time)
) %>% group_by(
Time, Species
) %>% summarise(
Richness = round(mean(Richness)),
.groups = "drop"
),
table(Species, Richness))
Richness
Species 0 1 2 3 4 5 6 7 8 9 10
1 80793 0 0 0 0 0 0 0 0 0 0
2 80623 0 0 0 0 0 0 1 0 0 169
3 80793 0 0 0 0 0 0 0 0 0 0
4 80047 0 0 0 0 0 0 0 2 0 744
5 79899 0 0 0 0 0 1 1 0 0 892
6 80793 0 0 0 0 0 0 0 0 0 0
7 80793 0 0 0 0 0 0 0 0 0 0
8 80793 0 0 0 0 0 0 0 0 0 0
9 80793 0 0 0 0 0 0 0 0 0 0
10 80793 0 0 0 0 0 0 0 0 0 0
11 80793 0 0 0 0 0 0 0 0 0 0
12 680 0 0 0 0 0 0 0 0 0 80113
13 80638 0 0 0 0 0 1 0 0 0 154
14 84 0 0 0 0 0 0 0 0 0 80709
15 80793 0 0 0 0 0 0 0 0 0 0
16 80793 0 0 0 0 0 0 0 0 0 0
17 80793 0 0 0 0 0 0 0 0 0 0
18 80793 0 0 0 0 0 0 0 0 0 0
19 80793 0 0 0 0 0 0 0 0 0 0
20 80793 0 0 0 0 0 0 0 0 0 0
21 80793 0 0 0 0 0 0 0 0 0 0
22 62 0 0 0 0 0 0 0 0 0 80731
23 80793 0 0 0 0 0 0 0 0 0 0
24 80545 0 0 0 0 1 0 0 0 0 247
25 80793 0 0 0 0 0 0 0 0 0 0
26 80793 0 0 0 0 0 0 0 0 0 0
27 80793 0 0 0 0 0 0 0 0 0 0
28 80793 0 0 0 0 0 0 0 0 0 0
29 80720 1 0 0 0 0 0 0 0 0 72
30 69383 0 1 1 1 0 0 2 1 2 11402
31 80793 0 0 0 0 0 0 0 0 0 0
32 500 0 0 0 0 0 0 0 0 1 80292
33 66478 0 0 1 3 4 1 1 0 2 14303
34 80793 0 0 0 0 0 0 0 0 0 0
35 80331 0 0 0 0 0 0 0 2 0 460
36 80793 0 0 0 0 0 0 0 0 0 0
37 46276 0 0 2 0 0 0 0 0 1 34514
38 80793 0 0 0 0 0 0 0 0 0 0
39 80793 0 0 0 0 0 0 0 0 0 0
40 80793 0 0 0 0 0 0 0 0 0 0
41 80793 0 0 0 0 0 0 0 0 0 0
42 80793 0 0 0 0 0 0 0 0 0 0
43 27626 0 0 1 0 2 0 1 1 0 53162
44 60221 0 0 0 0 3 1 0 0 3 20565
45 80793 0 0 0 0 0 0 0 0 0 0
46 80793 0 0 0 0 0 0 0 0 0 0
47 52181 0 0 1 0 1 1 0 0 0 28609
48 80793 0 0 0 0 0 0 0 0 0 0
49 80793 0 0 0 0 0 0 0 0 0 0
50 80793 0 0 0 0 0 0 0 0 0 0
51 80793 0 0 0 0 0 0 0 0 0 0
52 80793 0 0 0 0 0 0 0 0 0 0
53 80793 0 0 0 0 0 0 0 0 0 0
54 80793 0 0 0 0 0 0 0 0 0 0
55 78733 0 1 1 1 0 0 0 1 1 2055
56 80793 0 0 0 0 0 0 0 0 0 0
57 79019 0 0 0 0 0 0 0 0 0 1774
58 80793 0 0 0 0 0 0 0 0 0 0
59 80793 0 0 0 0 0 0 0 0 0 0
60 80793 0 0 0 0 0 0 0 0 0 0
61 80465 0 0 0 0 0 0 1 0 0 327
62 80793 0 0 0 0 0 0 0 0 0 0
63 80793 0 0 0 0 0 0 0 0 0 0
64 80793 0 0 0 0 0 0 0 0 0 0
65 80793 0 0 0 0 0 0 0 0 0 0
66 80793 0 0 0 0 0 0 0 0 0 0
67 80793 0 0 0 0 0 0 0 0 0 0
68 79018 0 1 0 1 0 0 1 2 4 1766
69 78129 0 1 0 0 2 0 2 1 0 2658
70 80793 0 0 0 0 0 0 0 0 0 0
71 1394 0 0 0 0 0 0 0 0 0 79399
72 292 0 0 1 0 0 0 0 0 0 80500
73 80793 0 0 0 0 0 0 0 0 0 0
74 80793 0 0 0 0 0 0 0 0 0 0
75 80793 0 0 0 0 0 0 0 0 0 0
76 78672 1 0 0 0 1 0 1 0 1 2117
77 79764 0 0 0 0 0 0 0 0 0 1029
78 79252 0 1 0 1 0 0 0 0 0 1539
79 80793 0 0 0 0 0 0 0 0 0 0
80 80793 0 0 0 0 0 0 0 0 0 0
81 80793 0 0 0 0 0 0 0 0 0 0
82 80793 0 0 0 0 0 0 0 0 0 0
83 80793 0 0 0 0 0 0 0 0 0 0
84 80793 0 0 0 0 0 0 0 0 0 0
85 80793 0 0 0 0 0 0 0 0 0 0
86 64857 0 0 1 1 3 0 1 1 1 15928
87 80793 0 0 0 0 0 0 0 0 0 0
88 80793 0 0 0 0 0 0 0 0 0 0
89 80793 0 0 0 0 0 0 0 0 0 0
90 80793 0 0 0 0 0 0 0 0 0 0
[ reached getOption("max.print") -- omitted 10 rows ]
AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.19145 0.11089 0.10844 0.07943 0.05993 0.04530
sum(summary(PCA)$importance[2, ])
[1] 0.99996
ggplot2::autoplot(PCA, loadings = TRUE,
data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::mutate(
dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCAPA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.21745 0.12196 0.08092 0.05942 0.05217 0.04471
ggplot2::autoplot(PCAPA, loadings = TRUE,
data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
sum(summary(PCA)$importance[2, ])
[1] 0.99996
ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)
Environments <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
retval <- data.frame(Environment = toString(i),
Time = time
)
retval <- cbind(retval, env)
colnames(retval) <- c("Environment", "Time",
paste0("Basal", 1:34),
paste0("Consumer", 35:100)
)
return(retval)
},
abund = result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
),
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)
EnvironmentDistance <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard")
EnvDistClust <- hclust(EnvironmentDistance)
plot(EnvDistClust, labels = FALSE)
EnvironmentDistancePA <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)
EnvDistClustPA <- hclust(EnvironmentDistancePA)
plot(EnvDistClustPA, labels = FALSE)